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ABSTRACT 

Context. Observations suggest that low-mass stars condense out of dense, relatively isolated, molecular cloud cores, with each core 
spawning a small-TV cluster of stars. 

Aims. Our aim is to identify the physical processes shaping the collapse and fragmentation of a 5.4 M G core, and to understand how 
these processes influence the mass distribution, kinematics, and binary statistics of the resulting stars. 

Methods. We perform SPH simulations of the collapse and fragmentation of cores having different initial levels of turbulence (a TURB = 
0.05, 0.10, 0.25). We use a new treatment of the energy equation that captures (i) excitation of the rotational and vibrational degrees 
of freedom of H2, dissociation of H2, ionisation of H and He, and (ii) the transport of cooling radiation against opacity due to both 
dust and gas (including the effects of dust sublimation, molecules, and H~ ions). We also perform comparison simulations using a 
standard barotropic equation of state. 

Results. We find that - when compared with the barotropic equation of state - our more realistic treatment of the energy equation 
results in more protostellar objects being formed, and a higher proportion of brown dwarfs; the multiplicity frequency is essentially 
unchanged, but the multiple systems tend to have shorter periods (by a factor ~ 3), higher eccentricities, and higher mass ratios. 
The reason for this is that small fragments are able to cool more effectively with the new treatment, as compared with the barotropic 
equation of state. We also note that in our simulations the process of fragmentation is often bimodal, in the following sense. The 
first protostar to form is usually, at the end, the most massive, i.e. the primary. However, frequently a disc-like structure subsequently 
forms round this primary, and then, once it has accumulated sufficient mass, quickly fragments to produce several secondaries. 
Conclusions. We believe that this delayed fragmentation of a disc-like structure is likely to be an important source of very low-mass 
stars in nature (both low-mass hydrogen-burning stars and brown dwarf stars); hence it may be fundamental to understanding the 
way in which the statistical properties of stars change - continuously but monotonically - with decreasing mass. However, in our 
simulations the individual cores probably produce too many stars, and hence too many single stars. We list the physical and numerical 
features that still need to be included in our simulations to make them more realistic; in particular, radiative and mechanical feedback, 
non-ideal magneto-hydrodynamic effects, and a more sophisticated implementation of sink particles. 
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1. Introduction 

The influence of thermodynamic effects on low-mass star for- 
mation in isolated, low-turbulence cores. There is direct and in- 
direct observational evidence to suggest that a significant frac- 
tion of low-mass stars form in small, relatively isolated, low- 
turbulence prestellar cores, with each core spawning only a few 
stars. This paper is concerned with exploring this mode of star 
formation, by means of numerical simulations, with a view to 
evaluating (a) how the statistical properties of the resulting pro- 
tostars depend on the initial conditions, and (b) what role is 
played by thermal, chemical, and radiative processes. We adopt 
an analytic approach. That is to say, we do not seek to implement 
all the deterministic physical effects at once, and hence we do 
not expect at this stage to reproduce all the observed features of 
real star formation. Rather, we seek to establish whether particu- 
lar physical effects, namely the thermal and radiative processes, 
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influence the outcome in a systematic way. Given the complex- 
ity of star formation and the limitations of numerical codes, this 
incremental approach seems a more fruitful one to pursue. 

The direct observational evidence for low-mass star forma- 
tion in isolated, low-turbulence cores. The direct evidence for 
this mode of star formation comes from a number of studies. 
First, Andre et al. (2007) have estimated the inter-core velocity 
dispersion in Ophiuchus (i.e. the dispersion in the bulk veloci- 
ties of cores relative to their neighbours). Andre et al. infer that 
the frequency of interactions between cores is so low that a typ- 
ical core is likely to have collapsed and fragmented, internally, 
before it interacts with a neighbouring core. Second, estimates 
of the level of turbulence in low-mass prestellar cores give sub- 
sonic values, i.e. cr TURB < cr THERM (e.g. Myers 1983; Myers et al. 
1991; Myers 1998; Andre et al. 2007). Indeed, Myers (1998) 
concludes that the decay of turbulence to subsonic levels may be 
a pre-requisite for the formation of low-mass protostars. 
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The indirect observational evidence for low-mass star for- 
mation in isolated, low-turbulence cores. Indirect evidence for 
this mode of star formation comes from the binary statistics 
of young low-mass stars, which show that a high fraction are 
born in binary or higher-order multiple systems. The binary frac- 
tion decreases with decreasing primary mass, and with increas- 
ing age, but for a 1 M primary it is still ~ 60% in the field 
(Duquennoy & Mayor 1991). Goodwin & Kroupa (2005) and 
Hubber & Whitworth (2005) have shown that this high multi- 
plicity requires newly-born stars to complete their early dynam- 
ical evolution in small sub-clusters containing just a few stars 
(i.e. Af SUBCLUSTER ~ 4 + 1 stars). This is because, in a small sub- 
cluster, A^-body interactions tend rather quickly to deliver a tight 
binary, usually comprising the two most massive stars, and to 
eject most of the remaining stars as singles (McDonald & Clarke 
1993; Goodwin & Kroupa 2005; Hubber & Whitworth 2005). 
Therefore, if AA SUBCLUSTER is increased, a higher proportion of stars 
are ejected as singles, and therefore the primordial binary frac- 
tion is reduced, in contradiction with the observations. 

The need for large ensembles of numerical simulations. One 
advantage of this mode of star formation, involving a low-mass 
core fragmenting in relative isolation to form a small number of 
protostars, is that it can be simulated with quite high numerical 
resolution. A related advantage is that many realisations can be 
computed. Given the chaotic nature of turbulent, self-gravitating 
gas dynamics, this is essential if robust statistical inferences are 
to be made. 

The effect of the initial level of turbulence in previous sim- 
ulations. The first comprehensive investigation of this mode of 
star formation was made by Goodwin et al. (2004a,b), who sim- 
ulated the evolution of a large ensemble of cores, all having the 
same mass, initial density profile and turbulent power spectrum. 
They considered three different levels of turbulence, character- 
ized by 

a TURB = 1 |sSL = 0.05, 0.10, 0.25, (1) 

where E TUSB and E ORAW are the initial turbulent and gravitational 
energies of the core. For each value of a TURB , they performed 
many different simulations, each with a different realisation of 
the turbulent velocity field, in order to obtain reasonable statis- 
tics. They found that increasing the level of turbulence increased 
both the total number of stars formed, and the proportion of 
brown dwarfs. 

The effect of the turbulent power spectrum in previous sim- 
ulations. In a subsequent paper, Goodwin et al. (2006) consid- 
ered different turbulent power spectra, of the form P k oc k~ x , 
and showed that increasing the exponent x resulted in more frag- 
ments. In effect, increased x means that there is more turbulent 
power on long wavelengths, and hence more large-scale frag- 
mentation resulting in separate protostars. This is in accordance 
with the findings of Klessen & Burkert (2001), who note that 
turbulence has two effects. Turbulence on very short (sub-Jeans) 
length-scales can be viewed as a source of extra pressure, supple- 
menting the thermal pressure. Conversely, turbulence on larger 
length-scales creates the structures which, if they become suffi- 
ciently massive and dense, are amplified by self-gravity to be- 
come prestellar cores. 

The limitations of using a barotropic equation of state. 
However, in their simulations Goodwin et al. (2004a,b, 2006) 
use a simple barotropic equation of state, which is designed 
to mimic the gross thermal behaviour of the gas at the centre 
of a collapsing, non-rotating 1 M protostar, as determined by 
Masunaga & Inutsuka (2000). This is not realistic. First, it does 



not take proper account of the thermal history of protostellar gas, 
which depends on environment, metallicity, mass and geometry. 
Second, a barotropic equation of state is unable to capture ther- 
mal inertia effects, for example, the complex system of interact- 
ing pressure waves which dominates the dynamics when the gas 
becomes adiabatic and the thermal timescale suddenly becomes 
longer than the dynamical timescale. This is when fragmenta- 
tion occurs (e.g. Boss et al. 2000), and thermal inertia effects are 
therefore critical. 

Our new treatment of the energy equation. We have therefore 
revisited the study of Goodwin et al. (2004a,b), but now using 
a significantly improved treatment of the energy equation due 
to Stamatellos et al. (2007a). This new treatment of the energy 
equation captures the critical thermal and radiative effects in pro- 
tostellar gas much more realistically than a barotropic equation 
of state. For full details of the new treatment, and of the tests 
that have been performed to establish its fidelity, the reader is 
refered to Stamatellos et al. (2007). With the new treatment of 
the energy equation, we find 

- that the efficiency of fragmentation is increased, in the sense 
that on average more stars are formed, and in particular more 
brown dwarfs; 

- that the mass distribution becomes bimodal, with a peak 
around 0.3 to 1 .0 M , and a subsidiary peak in the brown 
dwarf regime around 0.03 to 0.06 M (but note that we are 
only considering the fragmentation products of a single core 
mass - we do not expect this bimodality to translate to the 
overall stellar mass function, when a distribution of core 
masses is considered); 

- that fragmentation frequently involves the formation of 
a primary star (from material having low angular mo- 
mentum) and the subsequent accumulation of a massive 
disc around this primary, with the disc then fragmenting 
~20, 000 to ~ 100, 000 years later, to produce secondaries; 

- that brown dwarfs result from disc fragments which either 
are born in the outer disc where there is not much material 
to assimilate, or are quickly ejected from the disc, usually by 
mutual interactions, before they can assimilate much mass; 
disc fragments that form and remain in the inner disc usually 
assimilate sufficient mass to become hydrogen-burning stars; 

- that a significant fraction of the stars end up in binary sys- 
tems, or higher multiples, and these systems tend to have 
large mass ratios (i.e. q = M 2 /M 1 ~ 1); 

- that there is little evidence for competitive accretion - in fact, 
as noted above, there is a rather egalitarian process at work 
tending to produce binary systems with components of com- 
parable mass. 

The plan of the paper. In Section 2, we describe the initial 
conditions, the constitutive physics, and the numerical method. 
In Section 3, we present the results of our simulations and the 
statistical distributions derived from them, and relate these dis- 
tributions to the underlying physics. In Section 4 we summarise 
our conclusions. 

2. Simulations 

2. 1 . Initial conditions 

We use the same global initial conditions as Goodwin et al. 
(2004a,b, 2006). These initial conditions are intended to fit the 
observed properties of prestellar cores like L1544 (e.g. Ward- 
Thompson et al. 1994; Andre et al. 1996; Ward-Thompson et al. 
1999; Andre et al. 2000; Alves et al. 2001; Kirk et al. 2005). 
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Fig. 1. Two simulations of the collapse and fragmentation of a 5.4 M core having the same initial turbulent velocity field. In the 
first simulation (top row) the core is evolved using the barotropic equation of state. In the second simulation (bottom row) the 
core is evolved using the new treatment of the energy equation. Each frame shows the logarithm of the column density through the 
computational domain. For each simulation we show an image at, reading from left to right, t — 67 kyr , 73 kyr , 74 kyr , and 75 kyr . 
We see that with the new treatment of the energy equation, the disc is more unstable against fragmentation. 



Density profile and mass. The density structure of a low- 
mass prestellar core usually appears to consist of a central kernel 
having approximately uniform density, surrounded by an outer 
envelope in which the density falls off radially with exponent 
2 < \d£og(p) I d(og(r)\ < 5 . We therefore adopt a Plummer-like 
initial density profile, 



p{r) = 



(l + (r/fl KERNEL ) 2 ) 



(2) 



= 3 x 10 



18 gem 3 is the central density, and 



Here Pke 

^kernel = 5,000AU is the radius of the central kernel (within 
which the density is approximately uniform). The mass inside 
^kernel i s tnen ^kernel = 11 Mq. The outer envelope of the core 
extends out to R CORE = 50, 000 AU, so the total core mass is 
M C0RE = 5.4 M and the density at the boundary of the core is 
10 4 times lower than at the centre. 

Temperature. The gas is initially isothermal at T = 10 K, and 
hence the initial ratio of thermal to gravitational energy is 



0.3. 



(3) 



Turbulence. We also impose an initial divergence-free 
Gaussian random velocity field on each core. We set the power 
spectrum of this velocity field to be P^dk oc ]c A dk, since this 
mimics the scaling laws observed in molecular clouds (e.g. 
Larson 1981; Burkert & Bodenheimer 2000). This prescription 
for initialising the velocity field is normally referred to as tur- 
bulence (c.f. Bate et al. 2002a,b, 2003; Bonnell et al. 2003; 
Goodwin et al. 2004a,b, 2006), and we shall follow this con- 
vention. However, we emphasise that this is not self-consistent, 
fully developed turbulence (c.f. Offner et al. 2008), it is simply a 
device used to seed the core with substructure, which may then 
be dissipated or become amplified by self-gravity. In the simula- 
tions presented here, we consider quite low levels of turbulence, 
characterised by a TURB = 0.05, 0.10, 0.25 (where a TURB is defined 
in Eqn.[TJ. We note that these are typical values for the level of 



turbulence in observed low-mass cores, as collated in the cata- 
logue of Jijina, Myers & Adams (1999); they are much lower 
than the value a TVRB = 1 used by Bate et al. (2002a,b, 2003) and 
Bonnell et al. (2003) in their simulations of more massive cores. 

2.2. Constitutive physics 

Barotropic equation of state. For comparison purposes, we fol- 
low Goodwin et al. (2004a,b, 2006) in using a barotropic equa- 
tion of state, such that the isothermal sound speed, c s , is given 
by 



Clip) = 

P 



Po 



2/3 



(4) 



Here c Q = 0.2kms _I and p Q = 10~ 13 gcirT 3 . At low densities 
(p < p ) the gas is approximately isothermal; c Q = 0.2 km s -1 
corresponds to a mix of X — 0.70 molecular hydrogen and 
Y = 0.28 atomic helium at T = 10 K. At higher densities 
(p > p a ) the presumption is that the gas becomes opaque to 
its own cooling radiation, and heats up adiabatically. Until the 
temperature rises above ~ 200 K, the rotational degrees of free- 
dom of molecular hydrogen are not significantly excited, so the 
effective adiabatic exponent is y ~ 5/3 and the sound speed rises 
as ~p'^ 3 . This in turn means that the Jeans mass increases rather 
rapidly once the density exceeds p CRIT , roughly as M JEANS oc p'^ 2 . 
We note that the simulations of Bate et al. (2002a,b, 2003) and 
Bonnell et al. (2003) use a similar barotropic equation of state to 
us, but with y = 7/5 in the adiabatic regime; consequently their 
Jeans mass increases much more slowly with increasing density, 
M JEANS oc p 1 ^ 10 , giving a greatly extended window of opportunity 
for fragmentation at low masses. 

The need for a more realistic treatment of the energy equa- 
tion. Eqn. (0 is a reasonably good fit to the run of temperature 
with density at the centre of a collapsing, non-rotating 1 M Q pro- 
tostar, as obtained by the detailed computations of Masunaga & 
Inutsuka (2000). It also matches reasonably well the earlier re- 
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suits of Larson (1969) and Tohline (1982) for the thermodynam- 
ics of protostellar matter. However, even if one limits consid- 
eration to the collapse of a spherically symmetric, non-rotating 
1 M protostar, the run of temperature with density away from 
the centre shows quite a large variance, relative to the central val- 
ues. Moreover, as soon as the simulations involve condensations 
which have masses below 1 M and/or are non-spherical, these 
condensations will tend to become opaque, and to heat up, at a 
significantly higher density. This is because in general (a) the 
optical depths involved are lower, and (b) the rates of compres- 
sional heating are slower. A barotropic equation of state cannot 
capture these effects, nor can it handle situations where the ther- 
mal timescale becomes comparable with the dynamical one. Yet 
these are precisely the circumstances under which fragmentation 
occurs, as evidenced by the simulations of Boss et al. (2000). It is 
therefore of great interest to explore what happens when the en- 
ergy equation and associated radiative transport are treated more 
realistically. 

New treatment of the energy equation. Stamatellos et al. 
(2007a) have introduced a new treatment of the energy equa- 
tion and the associated radiative transport, which captures all the 
above effects, but incurs a very small computational overhead 
(~ 3%). This method incorporates the internal energy associ- 
ated with the rotational and vibrational degrees of freedom of 
H 2 , the dissociation of H 2 , and the ionisation of H°, He and 
He + . It also treats the transport of heating and cooling radiation, 
using an opacity which - in different regimes - is delivered by 
refractory dust cores with ice mantles, refractory dust cores with- 
out ice mantles, molecules, bound/free and free/free transitions 
and electron scattering. The background radiation field and the 
metallicity are user-defined free parameters. In the current ver- 
sion it is assumed that the dust abundance is proportional to the 
metallicity and that the dust properties do not evolve, but these 
assumptions can easily be relaxed. The method adapts efficiently 
to circumstances in which the thermal timescale is shorter than 
the dynamical one, and vice versa. The method is able to handle 
optically thin, intermediate and optically thick situations, and it 
has been extensively tested against detailed computations (like 
those of Masunaga & Inutsuka 2000) and analytic solutions (like 
those of Spiegel 1957). For full details of the method, and of 
the tests to which it has been subjected, the reader is referred to 
Stamatellos et al. (2007a). 

2.3. Numerical method 

The SPH code. The simulations are performed using the dragon 
smoothed particle hydrodynamics code (Goodwin et al. 2004a). 
This is a standard SPH code (e.g. Monaghan 1992), and has been 
extensively tested and optimised. It uses a second-order Runge- 
Kutta integration scheme, multiple-particle timesteps evalu- 
ated according to a Courant-Friedrich type condition, adaptive 
smoothing lengths and standard artificial viscosity (Monaghan 
1992). An octal spatial decomposition tree (Barnes & Hut 1986) 
is used to facilitate the collation of neighbour lists and the com- 
putation of gravitational accelerations. The smoothing lengths of 
the particles are adjusted so that each particle always has exactly 
N NB]B = 50 neighbours. In this way numerical diffusion is kept 
very low (Attwood et al. 2007). 

Sink particles. During SPH simulations of star formation, the 
timestep required to follow the motion of the gas decreases as the 
stellar condensations increase in density, and eventually all the 
computational resources are being used to follow the first one 
or two condensations to form. The rest of the simulation then 
ends up in a state of suspended animation, and so it is impos- 



sible to know whether any further stars would have formed, or 
what their properties might have been. In order to overcome this 
problem, we invoke sinks (Bate et al. 1995). Specifically, a sink 
is created if (a) the density of one of the SPH particles, i, ex- 
ceeds the threshold p SINK = 10" 11 gem -3 ; (b) particle ;' has neg- 
ative velocity divergence; and (c) particle i and its neighbours 
have net negative energy (i.e. are bound). The newly-formed 
sink incorporates all the neighbours of particle ;' within a dis- 
tance R smK = 5 AU. Any SPH particle which subsequently passes 
within a distance R slNK of a sink and is bound to that sink, is as- 
similated by the sink. We record and update the total mass, po- 
sition, velocity and angular momentum of each sink. The sinks 
created in a simulation are identified as stars. During post-run 
analysis, we can use their masses, positions and velocities to de- 
termine the mass distribution, kinematics and binary statistics of 
stars. 

Resolution. All the simulations reported here are performed 
with N TOT = 25, 000 SPH particles. Hence the mass resolution is 



Af NFTR M TnT 
M Mm ~ \f ~ °- 01M °' 

* * TOT 



(5) 



(Bate & Burkert 1997), and the formation of stars below this 
mass is inhibited (Whitworth 1998). Since sink particles with 
radius R smK = 5 AU are created above a threshold density, p SINK = 
10" 11 gem -3 , the minimum linear resolution is R smK = 5 AU. 
The discs which form in our simulations typically have radius 
R msc ~ 50 AU and half-thickness Z DISC ~ 5 AU, so they are not 
resolved in the vertical direction. However, the fragmentation of 
a disc is essentially a two-dimensional process, in the sense that 
the forces which drive the accumulation of matter into a proto- 
fragment (the forces which enter into the Toomre criterion) are 
in the plane of the disc. Therefore we do not believe that the 
limited resolution is critical to our conclusions (but see Nelson 
2006). 



2.4. Nomenclature 

Stars. We will refer to all sinks as stars, irrespective of their 
mass. The implication is that all objects which form by gravita- 
tional instability, on a dynamical timescale, are ultimately stars, 
irrespective of whether they have sufficient mass to burn hydro- 
gen or deuterium. We note that for gas with isothermal sound 
speed c s ~ 0.2 km s -1 the dynamical timescale (i.e. the freefall 
time for a marginally Jeans unstable fragment) is rather short, 



GM M 

60kyr 

10c 3 y I Mo 



(6) 



particularly for very low masses. 

Planets. In contrast, objects which form by core accretion 
of solid matter (and possibly the subsequent acquisition of a 
gaseous envelope), on a much longer timescale (typically more 
than 1 Myr), are planets. Such objects cannot form in the present 
simulations, because the physics of dust settling and aggregation 
is not included, and the duration of the simulations is too short. 

The Stellar Initial Mass Function. In the context of star 
formation it is entirely appropriate to discuss - as a single 
group - all objects which form by gravitational instability. It 
is appropriate because the Initial Mass Function (IMF) ap- 
pears to be continuous across the hydrogen-burning limit at 
~ 0.080 M Q (e.g. Chabrier 2003), and there is no evidence for 
its not also being continuous across the deuterium-burning limit 
at ~ 0.013 M Q (just very poor statistics). It would therefore seem 
both more sensible and more practical to refer to this distribution 
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of mass as the Stellar IMF, rather than the IMF for Stars, Brown 
Dwarfs and Planetary-Mass Objects that Form by Gravitational 
Instability. The really critical issue here is that, in the mass in- 
terval 0.001 to 0.010 Mo, there are probably both objects that 
should be called stars and objects that should be called planets; 
distinguishing them observationally could be rather hard. 

The continuity of statistical properties across the hydrogen- 
burning limit. It is also appropriate because all the other observa- 
tional evidence suggests that the statistical properties of young 
brown-dwarf stars and low-mass hydrogen-burning stars (i.e. 
kinematics, clustering properties, binary statistics, frequency 
and lifetimes of circumstellar discs, accretion rates and out- 
flow rates, etc.) form a continuum (e.g. Burgasser et al. 2007; 
Luhmann et al. 2007; Whitworth et al. 2007). 

Continuity of formation mechanisms across the hydrogen- 
burning limit. This continuity of statistical properties suggests 
that the bulk of a star's mass is committed to the star it ends up 
in before it discovers which nuclear fuels it can burn and which 
it cannot. In other words, it seems likely that the mechanisms 
which form hydrogen-burning stars of mass ~ 0.090 M also 
work to form brown-dwarf stars of mass ~ 0.070 M , and vice 
versa. 

Trends across the hydrogen-burning limit. That is not to say 
that there are not strong trends in the statistical properties of 
stars as the mass decreases to very low values. Indeed, one of 
the main conclusions of this paper is that as one crosses the hy- 
drogen burning limit towards lower masses, (a) the proportion of 
stars that are formed by disc fragmentation is steadily increasing, 
and (b) this is reflected in their statistical properties. 

2.5. Ensembles of simulations for different treatments of the 
thermodynamics and different initial levels of turbulence 

For each value of a TURB (= 0.05, 0.10, 0.25), we have performed 
an ensemble of 20 simulations using the barotropic equation of 
state (Eqn.|H, and an ensemble of 20 simulations using our new 
treatment of the energy equation (Stamatellos et al. 2007). This 
gives a total of 120 simulations. The 20 simulations in an en- 
semble — i.e. the set of 20 simulations all using the same treat- 
ment of the thermodynamics (barotropic equation of state or new 
treatment of the energy equation) and the same initial level of 
turbulence (same a TURB ) — are distinguished solely by having 
different realisations of the turbulent velocity field. This simply 
requires the use of a different random-number seed to generate 
the initial turbulent velocity field. 

Each simulation is evolved for a total of 300 kyr. For com- 
parison, the initial freefall time at the centre of the core is 
f FF - 40 kyr; and the first star usually forms after 50 to 70 kyr. 

3. Results 

Table Q] lists, for each simulation performed with the barotropic 
equation of state, the identifier (id); the initial level of turbulence 
(a TURB ), the total mass which ends up in stars (2{.M*}/M ), the 
total number of stars (Af*), and the total number of brown dwarfs 
(A/ BD ), at the end of the simulation (after 300 kyr); the types of 
multiple system that have formed (S = single, B = binary, T 
triple. Q = quadruple); and the masses of the individual stars 
(M*/M Q ; with a superscript indicating which ones are compo- 
nents of multiple systems). Table [2] lists the same information 
for the simulations performed using the new treatment of the en- 
ergy equation. To save space, we give here only the first twenty 
lines of each table (i.e. the results obtained with ff TURB = 0.05); 
the full tables are available online. 



3.1. Efficiency and timing of star formation 

The effect of varying the initial level of turbulence, when using a 
barotropic equation of state. Table[3]records - for each treatment 
of the thermodynamics and each initial level of turbulence - the 
number of different realisations (N REAL ), the efficiency (i.e. the 
mean fraction of the core mass converted into stars after 300 kyr, 
77 = Tj{M-t}/M C0RE ), and the mean number of stars formed from 
one core (AT*). With a barotropic equation of state, increasing 
the initial level of turbulence (a) decreases the efficiency of star 
formation, 77, and (b) increases the mean number of stars formed 
from a single core, A/*. The efficiency is reduced by increased 
turbulence, because the outer diffuse parts of the core are more 
vigorously dispersed, and therefore after 300 kyr they have not 
yet had time to fall back into the core and be incorporated into 
stars. An increased initial level of turbulence increases the to- 
tal number of stars formed from a single core, because it drives 
more vigorous local compression, and thereby creates more pro- 
tostellar seeds (i.e. more lumps which are sufficiently dense to 
be gravitationally unstable). It follows that the mean stellar mass 
decreases with increasing turbulence. 

The effect of varying the initial level of turbulence, when us- 
ing the new treatment of the energy equation. With the new treat- 
ment of the energy equation, these monotonic trends disappear. 
The efficiency, 77, and the mean number of stars formed, N+, are 
only weakly dependent on the initial level of turbulence. 

The effects of the new treatment of the energy equation. On 
the other hand, switching from the barotropic equation of state to 
the new treatment of the energy equation reduces the efficiency, 
j], somewhat, and significantly increases the mean number of 
stars formed from a single core, N+. There are two physical ef- 
fects at work here. 

Lower stellar masses. First, the new treatment of the energy 
equation promotes the condensation of very low-mass stars, by 
taking proper account of the thermal history and environment 
of the gas. With the new treatment of the energy equation, a 
small proto-fragment tends to be cooler, and thereby more in- 
clined to condense out. This is because the new treatment takes 
account of the fact that, being of lower mass (and probably 
also non-spherical), the column-density inhibiting the cooling 
of a low-mass proto-fragment is lower; and because it contracts 
more slowly, its heating rate is lower. Consequently its temper- 
ature is lower - at a given density - than the one prescribed by 
the barotropic equation of state, since the latter is based on the 
behaviour at the centre of a collapsing, spherical, non-rotating 
1 Mq protostar. This is illustrated in Fig. [T] which shows ex- 
actly the same initial conditions, evolved first with the barotropic 
equation of state (top frames), and then with the new treatment 
of the energy equation (bottom frames). It is evident that disc 
fragmentation is far more advanced in the simulation using the 
new treatment of the energy equation. 

More ejections and lower efficiency. Second, because in the 
simulations performed with the new treatment of the energy 
equation the number of stars formed is higher, but they are of 
lower mass, these stars are less effective at mopping up the resid- 
ual gas (thereby increasing their mass), and more likely to be 
ejected by dynamical interactions with other stars. Hence the 
amount of mass converted into stars is somewhat reduced, and 
the efficiency, r\, is lower. 

Disc fragmentation. In fact, there is a common pattern of 
star formation in many of these simulations, irrespective of the 
treatment of thermodynamics. The low angular momentum ma- 
terial in the core collapses quickly to form the first star (here- 
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Table 1. Results of the simulations performed using the barotropic equation of state with a TURB = 0.05, 0.10 and 0.25, at time 
t = 0.3 Myr. 



ID 


TURB 


Z(M*)/M 


A4 


^BD 


Mult 


MJM a 


N071 


0.05 


3.731 


4 





T 


1.297', 0.970, 0.749', 0.715' 


N072 


0.05 


3.867 


5 





T 


1.638', 1.204', 0.472', 0.285, 0.268 


N073 


0.05 


3.282 


5 





B 


1.216*, 1.111', 0.386, 0.186, 0.383 


N074 


0.05 


4.001 


7 


3 


Q 


1.174?, 1.106?, 0.811', 0.806", 0.049, 0.031, 0.024 


N075 


0.05 


3.307 


4 





T 


0.956', 0.911', 0.728', 0.712 


N076 


0.05 


3.915 


7 


2 


Q 


1.001*, 0.914", 0.694, 0.593", 0.583", 0.088, 0.042 


N077 


0.05 


3.646 


1 





s 


3.646 


N078 


0.05 


3.814 


6 





T 


0.998', 0.892', 0.854', 0.763, 0.199, 0.108 


N079 


0.05 


3.959 


3 





T 


2.319', 0.823', 0.817' 


N080 


0.05 


3.700 


1 





S 


3.700 


N081 


0.05 


3.690 


1 





s 


3.690 


N082 


0.05 


3.928 


5 


2 


T 


1.322', 1.286', 1.214', 0.070, 0.036 


N083 


0.05 


3.905 


2 





B 


2.459*, 1.446* 


N084 


0.05 


3.987 


2 





B 


3.056*, 0.931* 


N085 


0.05 


3.911 


2 





B 


2.151*, 1.760* 


N086 


0.05 


3.774 


6 


1 


Q 


1.037", 1.007", 0.741", 0.732", 0.219, 0.038 


N087 


0.05 


3.404 


1 





S 


3.404 


N088 


0.05 


3.874 


1 





s 


3.874 


N089 


0.05 


3.491 


5 





B 


1.005*, 0.934*, 0.695, 0.693, 0.164 


N090 


0.05 


3.778 


1 





S 


3.778 



Table 2. Results of the simulations using the new treatment of the energy equation with a TURB = 0.05, 0.10 and 0.25, at time 
/ - 0.3 Myr. 



ID 


^TURB 


Z{M*}/M Q 


A& 


^BD 


Mult 


MJM e 


T071 


0.05 


3.161 


8 


3 


B 


0.825*, 0.810*, 0.622, 0.494, 0.346, 0.029, 0.020, 0.015 


T072 


0.05 


2.212 


6 


1 


BT 


0.750*, 0.603*, 0.281', 0.279', 0.275', 0.024 


T073 


0.05 


3.200 


5 


1 


B 


0.877*, 0.870*, 0.696, 0.678, 0.079 


T074 


0.05 


3.561 


13 


6 


BT 


0.830*, 0.828*, 0.439', 0.432', 0.427', 0.328, 0.154, 0.039, 0.034, 0.023, 0.013, 0.007, 

0.007 


T075 


0.05 


3.252 


7 


2 


BB 


0.753* 1 , 0.748* 1 , 0.599* 2 , 0.553, 0.546* 2 , 0.032, 0.021 


T076 


0.05 


3.918 


9 


5 


Q 


1.118", 1.104", 0.722", 0.715", 0.079, 0.070, 0.056, 0.029, 0.025 


T077 


0.05 


3.884 


7 


3 


B 


1.338*, 1.149*, 0.702, 0.635, 0.032, 0.017, 0.011 


T078 


0.05 


3.559 


12 


3 


T 


0.679, 0.678', 0.505', 0.478', 0.435, 0.261, 0.191, 0.168, 0.087, 0.040, 0.023, 0.014 


T079 


0.05 


3.410 


10 


1 


T 


0.733', 0.722, 0.487', 0.481', 0.461, 0.186, 0.139, 0.095, 0.088, 0.018 


T080 


0.05 


3.894 


4 


1 


B 


1.480*, 1.237, 1.159*, 0.018 


T081 


0.05 


3.613 


5 


1 


B 


1.158*, 1.059*, 0.686, 0.645, 0.065 


T082 


0.05 


3.741 


7 


2 


T 


0.907, 0.884', 0.758', 0.754', 0.399, 0.026, 0.013 


T083 


0.05 


3.739 


5 


1 


B 


0.999, 0.983*, 0.912*, 0.815, 0.030 


T084 


0.05 


3.868 


8 


3 


B 


1.219*, 1.191*, 1.123, 0.139, 0.092, 0.051, 0.031, 0.022 


T085 


0.05 


3.749 


7 


2 


B 


0.882*, 0.850, 0.832, 0.804*, 0.332, 0.026, 0.023 


T086 


0.05 


3.898 


6 


1 


Q 


0.863", 0.850", 0.849", 0.840", 0.460, 0.036 


T087 


0.05 


3.434 


1 





s 


3.434 


T088 


0.05 


3.381 


8 


3 


BB 


1.435", 1.073*', 0.332* 2 , 0.320* 2 , 0.164, 0.021, 0.019, 0.017 


T089 


0.05 


3.259 


9 


2 


B 


0.645*, 0.617*, 0.612, 0.542, 0.342, 0.314, 0.146, 0.022, 0.019 


T090 


0.05 


3.811 


7 


2 


T 


2.100', 1.032', 0.248', 0.187, 0.122, 0.058, 0.034 



after the primary) on a timescale of 50 to 70kyr, i.e. a bit longer 
than the initial freefall time at the centre of the core. Then ma- 
terial with higher angular momentum forms a circumstellar disc 
around the primary. This circumprimary disc grows in mass (the 
rate of infall onto the disc is greater than the rate at which mass 
accretes from the inner disc onto the primary) until the disc be- 
comes Toomre unstable and fragments to form multiple secon- 
daries. The delay between the formation of the primary and frag- 
mentation of the circumprimary disc is typically between 20 and 
100 kyr. During this time the disc is accumulating mass. Once 



the disc becomes Toomre unstable it normally fragments to pro- 
duce between 3 and 5 stars, in the space of a few kyr. 

Accretion histories. This pattern of fragmentation is illus- 
trated on Fig. [2] where, for a selection of simulations, we plot 
stellar masses as a function of time. On most of these plots, 
we see the primary forming, then a delay whilst the circum- 
primary disc builds up, and finally - when the circumprimary 
disc becomes Toomre unstable - the formation of a clutch of 
secondaries. Some of these secondaries are quickly ejected, and 
therefore end up as brown dwarfs, but others remain in the disc 
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Table 3. Basic statistics for each treatment of the thermodynamics (barotropic equation of state or new treatment of the energy 
equation) and each value of the initial level of turbulence (or^g). 



THERMODYNAMICS 


"turb 




V 


N-,, 


S 


B 


T 


Q 


mf 


cp 


cf 


BAROTROPIC EQUATION OF STATE 


0.05 


20 


0.694 


3.45 


29 


5 


6 


3 


0.33 


0.58 


1.19 




0.10 


20 


0.658 


4.15 


30 


2 


7 


7 


0.35 


0.64 


1.57 




0.25 


20 


0.600 


5.05 


27 


2 


10 


10 


0.46 


0.73 


1.82 


NEW TREATMENT OF THE ENERGY EQUATION 


0.05 


20 


0.605 


7.20 


88 


15 


6 


2 


0.21 


0.39 


0.63 




0.10 


20 


0.629 


6.20 


67 


7 


9 


4 


0.23 


0.46 


0.94 




0.25 


20 


0.623 


9.05 


98 


12 


13 


5 


0.23 


0.46 


0.90 




Fig. 2. Stellar masses as a function of time, for a selection of simulations. Note (i) the delay between the formation of the primary 
and the formation of a clutch of secondaries (this is the time during which the circumprimary disc accumulates, until it becomes 
Toomre unstable); and (ii) the marked decline in the accretion rate onto the primary once the secondaries start to condense out. 



and accrete sufficient mass to become hydrogen-burning stars. 
Occasionally some even grow bigger than the primary. 

The growth-time and fragmentation-time for the disc. In Fig. 
[3] we show the distributions of f, (the time of formation of the 
first star); f 2 - f, (the delay between the formation of the first 
and second stars); and f 3 - f 2 (the delay between the formation of 
the second and third stars), f, is the time it takes the low an- 
gular momentum material to assemble into the first stars and 
should be compared with the freefall time at the centre of the 
core (~ 40kyr); it has mean fi h = 61 kyr and standard deviation 
cr n = 6.2 kyr. (f 2 - f, ) is the time it takes to assemble a Toomre- 
unstable disc around the primary; it has mean fifo-to = 20 kyr 
and standard deviation cr^-to - 20 kyr, so there is quite a large 
range in the times required to assemble a Toomre unstable disc. 
(f 3 - f 2 ) is the time delay between the formation of the first two 
disc fragments; it has mean H(t 3 -t 2 ) - 1-6 kyr and standard devi- 
ation <T(f 3 -t 2 ) = 1.4 kyr. This short mean (?3 - t-i) reflects the fact 
that when the disc becomes unstable it becomes unstable over 
quite a large area, and therefore it tends to spawn several stars in 
quick succession. 

The effect of Toomre instability on the growth of the pri- 
mary. Another common feature of the accretion histories is that, 
when the circumprimary disc becomes Toomre unstable and 
fragments, the accretion rate onto the primary declines rapidly. 
Material which up until this juncture had been spiraling inwards 



and onto the primary star is now being used to create secondaries 
in the disc. 

Birth order. Fig.[4]shows the final mass of every star plotted 
against its formation time. Although the more massive stars tend 
to form earlier, the correlation is fairly weak. In all cases there 
is a delay before any brown dwarfs form. This is because most 
of the brown dwarfs, and also some of the low-mass hydrogen- 
burning stars, form in discs around more massive stars, and these 
discs take time to accumulate. 



3.2. The mass distribution of protostars 

Lower protostellar masses with the new treatment of the energy 
equation. Material which is parked in a circumprimary disc has 
time to lose entropy - to an extent that material which is com- 
pressed impulsively by turbulence does not. Consequently the 
masses of disc fragments are low, as predicted by Whitworth 
& Stamatellos (2006), and demonstrated by detailed numeri- 
cal simulations in Stamatellos et al. (2007b) and Stamatellos & 
Whitworth (2008a,b). However, this effect can only be captured 
with the new treatment of the energy equation, since this treat- 
ment takes account of the slow rate of compressional heating for 
matter parked in the disc, and the relatively low local column- 
densities through which its cooling radiation has to diffuse. In 
contrast, the barotropic equation of state presumes that the mat- 
ter is part of a spherical 1 M protostar, which by virtue of col- 
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Fig. 3. (a,d) The delay, f, , between the start of the simulation and the formation of the first star (the primary). (b,e) The delay, t^ - f, , 
between the formation of the first and second stars. (c,f) The delay, f 3 - t 2 , between the formation of the second and third stars. 
The top row (a,b,c) is for the simulations performed using the barotropic equation of state, and the bottom row (d,e,f) is for the 
simulations performed using the new treatment of the energy equation. 
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Fig. 4. The final mass of each star (at 300 kyr) against its formation time. The top row gives the results obtained with the barotropic 
equation of state for (a) a XURB = 0.05, open circles; (b) a mm - 0.10, open triangles; and (c) a XURB = 0.25, open stars. The lower row 



gives the results obtained with the new treatment of the energy equation for (d) a T 
triangles; and (f) a TURB = 0.25, filled stars. 



0.05, filled circles; (e) a 



0.10, filled 



lapsing more rapidly is heated more vigorously, and has to cool 
through a larger column-density; therefore, at a given density, it 
is hotter and fragments less readily (i.e. into more massive frag- 
ments, if at all). 



Mass distributions. The lower masses and greater numbers of 
stars formed with the new treatment of the energy equation pre- 
disposes the stars to mutual dynamical interactions which eject 
many of them before they have time to grow much by accretion. 
Fig- H shows the mass distributions obtained with the different 



R. E. Attwood, S. P. Goodwin, D. Stamatellos and A. P. Whitworth: Simulating star formation in molecular cloud cores 9 




ag[M] 
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Fig. 5. Normalised stellar mass distributions. The top row gives the mass distributions obtained with the barotropic equation of state 
for (a) a TURB = 0.05, (b) a Tl]RB = 0.10, and (c) a rum - 0.25. The lower row gives the mass distributions obtained with the new 
treatment of the energy equation for (d) a TURB = 0.05, (e) a TORB = 0.10, and (f) - 0.25. The black lines are histograms of the 
raw data, obtained using 15 equal logarithmic bins in the interval — 2 < log I0 (M*/M ) < + 1 , and the red lines are obtained by 
smoothing each protostellar mass with a Gaussian whose width is proportional to the separation between neighbouring masses (see 
Appendix A for details). 



combinations of thermodynamic treatment and initial level of 
turbulence. The black line shows the histogram obtained by dis- 
tributing the final stellar masses into 15 logarithmic bins which 
are equally spaced in the interval 



•2 < lot 



10 



A4 



< 1 



Alog 10 (M*) = 0.2. 



The red line shows the mass distribution obtained when each 
stellar mass is smoothed using a Gaussian smoothing kernel with 
adaptive smoothing lengths dictated by the separation between 
masses (see Appendix A for details). Both the histogram, and 
the smoothed distribution, are normalised, in the sense that 



f 

*J M=i 



dlog l0 (M) = 1. 



(7) 



From Fig. we see that the initial level of turbulence has little 
influence on the form of the mass distribution. 

Bimodal mass distributions. However, switching from the 
barotropic equation of state to the new treatment of the energy 
equation not only increases the proportion of brown-dwarf stars 
formed, but actually produces a bimodal mass distribution. The 
larger mode comprises hydrogen-burning stars with masses con- 
centrated in the range 0.3 to 1 .0 M , whilst the smaller mode 
comprises brown dwarf stars with masses concentrated in the 
range 0.02 to 0.06 M . This smaller mode represents very low- 
mass stars formed by disc fragmentation (due to the enhanced 
cooling which low-mass fragments enjoy with the new treatment 
of the energy equation) and then ejected by mutual interactions 
(before they can grow much by accretion). 



3.3. Multiplicity 

Measures of multiplicity. To discuss the statistics of multiplicity, 
we adopt the conventions proposed by Reipurth & Zinnecker 
(1993). We define systems to include single stars, and mul- 
tiple systems to include only systems that contain more than 
one star. If S is the number of single stars, B the number of 
binary systems, T the number of triple systems, Q the num- 
ber of quadruple systems, etc., then the total number of sys- 
tems is (S + B + T + Q + ...), the total number of multiple 
systems is (B + T + Q + ...), and the total number of stars is 
(5 + 2B + 3T + 4Q + ...). We can then compute the multiplicity 
frequency, mf, which measures the fraction of systems which are 
multiple, i.e. 



mf 



(B + T + Q + ...) 
(S +B + T + Q + ...)' 



(8) 



the companion probability, cp, which measures the fraction of 
stars which are in multiple systems, i.e. 



cp = 



(2B + 3T +4Q + ...) 
(S +2B + 3T + 4Q + ...) ' 



(9) 



and, from Goodwin et al. (2004b), the companion frequency, cf, 
which measures the mean number of companions which a star 
has (irrespective of whether it is a primary), i.e. 



cf = 



(2B + 6T + 12Q + ...) 
(S + 2B + 3T + 4Q + ...) 



(10) 



Multiplicity statistics. In Table [3] we record — for each en- 
semble of 20 simulations, representing a particular combination 
of thermodynamic treatment and initial level of turbulence — 
the total numbers of singles (S), binaries (B), triples {T) and 
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Fig. 6. For each multiple system we plot the number of stars 
formed in that simulation, A4, against the periods, P. (a) Results 
obtained using the barotropic equation of state; here open cir- 
cles represent aturb = 0.05, open triangles represent ff tm b = 0.10, 
and open stars represent a tur b = 0.25. (b) Results obtained us- 
ing the new treatment of the energy equation; here filled circles 
represent a tur b = 0.05, filled triangles represent arturb = 0.10, and 
filled stars represent a t urb = 0.25. 



quadruples (Q) formed in all simulations; and the mean mul- 
tiplicity frequency (mf), the mean companion probability (cp), 
and the mean companion frequency (cf). 

Computing periods. Fig. [6] shows, for each simulation, the 
number of stars formed in that simulation plotted against the 
periods of all the multiple systems identified at the end of the 
simulation. These periods are derived on the assumption that all 
multiple systems are hierarchical, which is not always true. Thus 
the two periods for a triple system are extracted by finding the 
period for the pair with the greatest specific binding energy, then 
treating this pair as a single star and finding the period of its orbit 
relative to the third star. This is appropriate for stable hierarchi- 
cal systems, but of limited value for unstable non-hierarchical 
systems. 

Subsequent evolution of periods. We should therefore expect 
some subsequent evolution in these distributions, with mutual 
interactions tending to lead to close systems becoming more 
tightly bound (occasionally with exchange of components) and 
wide systems being disolved. Eventually there will also be inter- 
actions with stars formed in neighbouring cores. These interac- 
tions will further disrupt the wider systems but have little effect 



on the closer systems. However, our simulations are not contin- 
ued long enough for such interactions to be important. 

Period distributions. The period distribution obtained using 
the barotropic equation of state has a mean /Ui ogl0 (p/ yl -) - 2.2 and a 
standard deviation <T logio( p/ yl ) 1 .0; periods range from ~10 4 yr 
down to ~ lOyr, and very few simulations form more than 7 
stars. With the new treatment of the energy equation, the periods 
are on average shorter (by about a factor of 3), with a mean of 
p\og m (Piyx) — 1.7, and a standard deviation o"iog 10 (p/yr) - 10; peri- 
ods range from ~ 10 4 yr down to ~ 3 yr, and many simulations 
form more than 7 stars. 

Effect of the new treatment of the energy equation on the 
period distribution. The reason why the new treatment of the en- 
ergy equation results in shorter-period multiples is that it allows 
the gas - in particular, the gas in smaller proto-fragments - to 
stay cooler to higher densities. Consequently the Jeans length, 
and hence the separations between neighbouring stars, tend to 
be smaller. 

Effect of the initial level of turbulence on the period distribu- 
tion. There is no obvious dependence of the period distribution 
on the level of turbulence, although this must be set against the 
poor statistics (between 26 and 53 periods for each combination 
of thermodynamics and initial level of turbulence). 

Unresolved orbits. We should also caution that the low- 
period systems are poorly resolved, in the sense that at periastron 
the stars are closer together than ^ SINK , and therefore their grav- 
itational interaction is softenned. This means that they should 
probably be somewhat more tightly bound. We have checked 
the formation of the individual stars in some of these close sys- 
tems, and established that in each case the two constituent stars 
(i.e. sinks) were initially created from well-defined and separate 
Jeans-unstable density peaks. 

Distribution of eccentricities. Fig.|7]shows orbital eccentric- 
ities (<?) plotted against periods (P), at the end of the simulations. 
The eccentricities are not strongly correlated with period, nor - 
modullo the poor statistics (see above) - do they appear to be 
correlated with the initial level of turbulence. However, there is 
a noticeable difference between the distributions obtained with 
the two different treatments of the thermodynamics. Using the 
barotropic equation of state, the distribution is concentrated to- 
wards high eccentricities, but there is still a substantial fraction, 
~ 25%, of systems having approximately circular orbits, e S= 0.2. 
Using the new treatment of the energy equation, the distribution 
of eccentricities is more strongly skewed towards high values, 
and less than 6% have e ^ 0.2. 

Effect of the treatment of thermodynamics on the distribu- 
tion of eccentricities. The barotropic equation of state facilitates 
the formation of low-eccentricity binaries by making it harder 
for circumbinary discs to fragment. At a given density the gas 
is hotter. Consequently, quite massive but relatively warm cir- 
cumbinary discs resist further fragmentation and instead act to 
dampen orbital eccentricities by accreting slowly onto the ex- 
isting binary components. In contrast, when the new treatment 
of the energy equation is used, massive circumbinary discs are 
relatively cool, so they fragment, and interactions between these 
additional fragments and the original components of the binary 
act to amplify the orbital eccentricities. 

Distribution of mass ratios. Fig.[8]shows the distributions of 
mass ratio, q = MJM X , at the end of the simulations. The dis- 
tributions are strongly skewed towards q ~ 1, i.e. nearly equal 
component masses. The mass ratios do not appear to be cor- 
related with the initial level of turbulence, a TURB , but again the 
statistics are poor. Mass ratios are correlated with orbital peri- 
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Fig. 7. Orbital eccentricities, e, plotted against periods, P, for 
multiple protostars: (a) Results obtained using the barotropic 
equation of state; here open circles represent a t urb = 0.05, 
open triangles represent or tur b = 0.10, and open stars represent 
o'turb = 0.25. (b) Results obtained using the new treatment of 
the energy equation; here filled circles represent a lm \, = 0.05, 
filled triangles represent a tur b = 0.10, and filled stars represent 
a t urb = 0.25. 



ods, in the sense that shorter-period systems tend to have higher 
mass-ratios, which is comparable with observations (e.g. Mazeh 
et al.1992). Since simulations conducted with the new treatment 
of the energy equation tend to produce multiples with shorter 
periods, they also tend to produce multiples with higher mass 
ratios. 

Mass equalisation in close binary systems. A mechanism 
which drives mass ratios towards unity in simulations of star 
formation was first described by Chapman et al. (1992), and has 
subsequently been noted by Burkert & Bodenheimer (1993) and 
by Bate & Bonnell (1997) (but see Ochi et al. 2005 for a different 
view, and Clarke 2007 for a rebuttal of Ochi et al.). If a binary 
system continues to grow by accretion, the specific angular mo- 
mentum of the infalling material (relative to the centre of mass 
of the binary system) tends to increase with time. Consequently 
the component with lower mass (the secondary, M 2 ), which nec- 
essarily is on a more extended orbit, is better disposed to assim- 
ilate this material with high angular momentum, and therefore 
it grows in mass until it is comparable with the primary (M,). 
This is the mechanism that appears to be operating here. It is 
less effective in wide binary systems, because the components 
of a wide system tend to accrete from separate reservoirs. 
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Fig. 8. The distribution of mass ratios, q, for multiple protostars: 
(a) using the barotropic equation of state; (b) using the new treat- 
ment of the energy equation. 



Competitive accretion and ejection. We note that there is lit- 
tle evidence for competitive accretion in our simulations. The 
first star to form (the primary) is more often than not the most 
massive at the end. However, stars forming later in the simu- 
lation frequently acquire comparable, and occasionally greater, 
masses than the primary. The material which ends up in these 
stars is normally rather coherently located. For example, once 
the circumprimary disc has formed, the material destined to form 
a particular secondary star accumulates in a particular range of 
radii, and sits there until it is mopped up by the growing sec- 
ondary star. Ejection does play a role in separating some stars 
from the reservoir of material they might otherwise have ac- 
creted, and thereby creating very low-mass stars. However, the 
material which accretes onto a star was in general present at the 
star's inception; its self-gravity contributed to the condensation 
which triggered the formation of a sink by pushing the density 
above p SINK . 

3.4. Missing physics 

The switch from the standard barotropic equation of state to our 
new more realistic treatment of the energy equation produces 
significant changes in the statistical properties of the stars re- 
sulting from the collapse and fragmentation of an isolated, low- 
turbulence, 5.4 M core. However, there are several important 
physical effects missing from our simulations. In particular, there 
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is no feedback from the stars, there are no (non-ideal) MHD ef- 
fects, and the use of sink particles raises some concerns. 
Feedback. Feedback from stars can take several forms. 

(i) The radiation from the stars will heat the surrounding dust 
and gas. Krumholz, Klein & McKee (2007) have recently simu- 
lated the collapse and fragmentation of more massive cores (100 
and 200 M ) with much higher initial levels of turbulence than 
those invoked in our simulations. Their treatment of the ther- 
modynamics takes account not only of the energy equation and 
the transport of cooling radiation, but also of radiative feedback 
from the forming stars in the core. They find that their cores only 
spawn a small number of stars. This is because the primary pro- 
tostar, which forms early on from material with relatively low 
angular momentum, has a high luminosity, and therefore sta- 
bilises the inner parts of its circumstellar accretion disc, by heat- 
ing them up. Fragmentation is only possible in the outer more 
diffuse parts of the disc. We expect a similar - but more mod- 
est - effect in low-mass, low-turbulence cores, more modest be- 
cause the primary luminosity will be much smaller. Nonetheless, 
it is likely that, even in the low-mass regime, the luminosity of 
the primary star is sufficient to inhibit fragmentation of the in- 
ner disc. The analytic work of Whitworth & Stamatellos (2006) 
predicts that a disc around a Sun-like primary is unlikely to 
fragment inside ~ 100 AU, and this is confirmed by the simula- 
tions of Stamatellos et al. (2007b) and Stamatellos & Whitworth 
(2008a,b). Consequently, the primary will end up more massive 
(by accreting the matter which is unable to fragment); the cir- 
cumprimary disc will take longer to grow to Toomre instability; 
and the secondaries which then condense out of it will be smaller 
in number, and at larger radii. 

(ii) Mechanical feedback, in the form of bipolar outflows 
will punch holes in the core. Preliminary exploration of this 
phenomenon (Stamatellos et al. 2005) suggests that it does not 
greatly change the efficiency of star formation, but it does slow 
it down (i.e. the delay between the formation of the primary and 
the formation of the secondaries is longer). This needs to be ex- 
plored further. 

(iii) Ionising radiation and winds from massive stars produce 
more violent feedback. We have recently developed the numeri- 
cal machinery to explore this (Bisbas et al., in preparation), but 
it is not part of the star-formation mode with which we are con- 
cerned here. 

MHD. Non-ideal MHD effects are likely to play an impor- 
tant role, and we have developed a code to treat them (Hosking 
& Whitworth 2004). However, it is a rather crude and ineffi- 
cient code, and further work is ongoing to improve it to the 
stage where it can be used to perform a large ensemble of sim- 
ulations. Price & Bate (2008) have simulated the collapse and 
fragmentation of a massive magnetised turbulent core, using an 
ideal MHD code, with a barotropic equation of state. They find 
that the magnetic field reduces both the efficiency of star forma- 
tion (i.e. the fraction of the core mass which ends up in stars) 
and the production of brown dwarfs. In an earlier paper (Price 
& Bate 2007), using more idealised initial conditions (a spheri- 
cal uniform-density cloud with an imposed m — 2 perturbation), 
they have shown that a magnetic field can also inhibit disc frag- 
mentation, by slowing the rate of disc growth and accelerating 
the rate at which angular momentum is redistributed. 

Sinks. Finally, we note that the use of sinks may com- 
promise our results, in ways which are hard to quantify (e.g. 
Commercon et al. 2008). First, the use of sinks means that all 
processes on scales below ~ 2R sm¥ _ = 10 AU, and at densities 
above p SINK = lO^gctrT 3 , are at best not properly resolved 
(e.g. orbits), and at worst excised completely (e.g. the second 



collapse when molecular hydrogen dissociates). Second, the cre- 
ation of sinks promotes TV-body interactions, and hence ejections 
of stars, whilst suppressing dissipative interactions between, and 
mergers of, stars. One can postpone the creation of sinks until 
very high densities are reached. For example, Stamatellos et al. 
(2007b) use p SINK = 10~ 2 g cirT 3 . However, this is very expensive 
computationally. This is an area in which a more sophisticated 
algorithm is needed. 

3.5. Comparison with observation 

Since we only treat a single core mass, with a single initial ra- 
dius and a single initial density profile, and since - as discussed 
in the preceding section - there are several potentially critical 
physical effects which are not included in our simulations, we 
do not expect to reproduce all the observed statistical properties 
of real stars. Nonetheless, it is appropriate to rehearse the vari- 
ous counts on which the properties of stars formed in our simu- 
lations conform to, or diverge from, reality; and to speculate on 
the reasons. 

Mean number of stars per core, N*. Our simulations form 
too many stars per core, and furthermore this over-production 
of stars is significantly exacerbated by the switch from the 
barotropic equation of state to the new, more realistic treatment 
of the energy equation. This is because the new treatment al- 
lows circumstellar discs, and low-mass fragments thereof, to 
stay cool to higher densities than the barotropic equation of state 
(which treats all gas as if it were at the centre of a collapsing, 
non-rotating 1 M protostar). The analytic results of Rafikov 
(2005), Matzner & Levin (2005), Kratter & Matzner (2006) 
and Whitworth & Stamatellos (2006), and the numerical simula- 
tions of Krumholz, Klein and McKee (2007), Stamatellos et al. 
(2007b) and Stamatellos & Whitworth (2008a,b) all suggest that 
the inclusion of radiative feedback reduces the number of stars 
formed, essentially by heating the inner disc, and thus suppress- 
ing Toomre instability by increasing the cooling time (Toomre 
1964; Gammie 2001). Whitworth & Stamatellos (2006) show 
that, given a solar-mass primary star at the centre of the disc, it 
can only fragment at large radii, R ^ 100 AU. The inclusion of 
mechanical feedback (Stamatellos et al. 2005) and/or a magnetic 
field (Price & Bate 2008) is also likely to reduce the number of 
stars formed, and in particular the number of brown dwarfs, by 
reducing the rate of accretion onto the primary and its circumpri- 
mary disc. Indeed, these effects are probably essential to reduce 
the efficiency of star formation in low-mass cores to the levels in- 
fered from observation. These levels are typically ~ 30 % (e.g. 
Nutter & Ward-Thompson 2007; Simpson et al. 2008). 

Mass distribution. The overall mass distribution produced by 
a single core, as a fraction of the core's total mass, is not con- 
strained by observation; if it were, we would know how to map 
the observed core mass function into the stellar initial mass func- 
tion. One interesting feature of our results is the suggestion that, 
amongst the stars spawned by a single core, there might be a bi- 
modal distribution of masses, comprising primary stars formed 
relatively early on, and secondary stars of much lower mass 
formed somewhat later by disc fragmentation. We also note that 
the mass of core that we are simulating (~ 5.4M ) is rather larger 
than the average isolated core (e.g. Alves et al. 2007; Nutter & 
Ward-Thompson 2007). This further complicates any attempts 
to map the results of these simulations onto the observed distri- 
butions which are the sum of a variety of core masses, mostly 
somewhat smaller than our core. 

Multiplicity. The multiplicity frequency of the stars formed 
in our simulations, mf ~ 0.2, is too low, especially for the 
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higher-mass stars; for the brown dwarfs and very low-mass 
hydrogen-burning stars (those with M+ < 0.1 M Q ) mf ~ 0.2 is 
actually in the middle of the range inferred from the limited ob- 
servations available (Burgasser et al. 2007; Luhman et al. 2007; 
Joergens (2008). The multiplicity frequency is expected to rise if 
the inclusion of extra physics reduces the number of stars formed 
from a single core. If this reduction is attributable to the suppres- 
sion of fragmentation in the inner parts of circumprimary discs, 
then the simulations of Stamatellos et al. (2008b) suggest that it 
will increase the multiplicity frequency of the higher mass stars 
(M^ ~ M ), and have little effect on the multiplicity frequency 
of the very low-mass stars (M* <0.1 M ); the simulations would 
then accord better with the observed distribution of multiplic- 
ity frequency, which appears to be a monotonically decreasing 
function of primary mass (e.g. Joergens 2008). 

Binary periods. The periods, P, of the binary systems formed 
in our simulations fall in the range 3 £ [P/yr] £ 10 4 . Systems 
with shorter periods cannot be resolved, because the gravita- 
tional fields of sink particles are softened at distances closer than 
R = 5 AU. Systems with longer periods must either form in more 
extended cores than the one we have modelled here, or they must 
result from interactions between stars formed in separate cores. 
An encouraging feature of the multiple systems formed in our 
simulations is the fact that most of the very low-mass systems 
(M 1 < 0.1 M Q ) have periods in the range 10 to lOOyr, in good 
agreement with the separations of observed very low-mass sys- 
tems (e.g. Joergens 2008). 

Mass ratios and eccentricities. The mass ratios of the mul- 
tiple systems formed in our simulations are concentrated to- 
wards high values. This again accords with what is observed 
for very low-mass systems (e.g. Burgasser et al. 2007; their Fig. 
5a), but contrasts with the flatter distribution observed in higher- 
mass systems. The multiple systems formed in our simulations 
are also skewed towards much more eccentric orbits than ob- 
served systems. However, the eccentricity distribution at birth 
is almost impossible to compare to the distribution in older sys- 
tems. Firstly, close systems will be circularised by tidal and other 
dissipative forces. Secondly, wider systems will be subject to 
encounters which will rapidly change the birth eccentricity dis- 
tribution beyond recognition (Parker et al., in preparation). Our 
simulations do not address these possibilities. 

3.6. Convergence 

We have repeated one of our simulations with 50,000 and 80,000 
sph particles, to check whether our simulations are converged, 
in a statistical sense (i.e. whether the statistical distributions of 
stellar parameters does not depend on the number of sph parti- 
cles used). For this purpose we have chosen simulation T011, 
which has an initial level of turbulence a TURB = 0.10, and uses 
the new treatment of the energy equation; the mean number of 
stars formed with this combination is N+ = 6.2 (see Table [3]). 
In the original T011 simulation, with just 25,000 sph particles, 
6 stars are formed. With 50,000 sph particles 6 stars are again 
formed. With 80,000 sph particles 8 stars are formed. We stress 
that in this context convergence can only be discussed in a sta- 
tistical sense. This is because, with the low initial levels of tur- 
bulence we are using, the gravitational fragmentation that en- 
sues is seeded from two sources. There are the small density 
enhancements created by subsonic converging flows due to the 
initial imposed turbulent velocity field; these are reproducable 
when using different particle numbers. However, there is also 
particle noise; this is not reproducable when using different par- 



ticle numbers. Therefore convergence can only be tested fully 
by repeating the whole ensemble of simulations with higher par- 
ticle numbers, and this is not feasible with the computational 
resources at our disposal. We are currently preparing a paper 
which demonstrates convergence in a simulation of gravitational 
fragmentation by using very carefully relaxed initial conditions; 
the imposed perturbation (which is reproducable) is then able to 
dominate particle noise in seeding gravitational fragmentation. 
These simulations exhibit excellent convergence, as do the sim- 
ulations of Jeans instability presented by Hubber et al. (2006). 
We are therefore confident that our code is capturing gravita- 
tional fragmentation faithfully. 

4. Conclusions 

We have performed a large ensemble of SPH simulations of 
the collapse and fragmentation of an isolated, turbulent 5.4 M 
core, with a view to establishing how the statistical properties 
of the resulting stars are influenced by (i) different initial levels 
of turbulence, and (ii) different treatments of the thermodynam- 
ics. We consider three initial levels of turbulence, characterised 
by a TURB = 0.05, 0.10 and 0.25. We treat the thermodynamics 
firstly with a standard barotropic equation of state, and secondly 
with a new treatment of the energy equation which captures all 
the important energy modes of the gas and takes account of ra- 
diation transport and variations in the opacity. 

Increasing the initial level of turbulence tends to reduce the 
efficiency of star formation, r\ (i.e. the fraction of the core mass 
which is converted into stars after 300 kyr), and to increase the 
number of stars formed by a single core, A4 , but the effect is very 
small, and all the other statistical properties of the stars formed 
are essentially independent of a TURB . 

A common pattern is observed in which the low-angular- 
momentum material in the core collapses to form the first star 
(the primary) after 50 to 70 kyr (i.e. just a bit longer than the ini- 
tial freefall time at the centre of the core), and then a massive disc 
builds up around the primary. As soon as this circumprimary disc 
becomes Toomre unstable (which may take from 20 to 100 kyr), 
it rapidly breaks up into a clutch of secondary stars. Those sec- 
ondaries that are quickly ejected from the disc normally end up 
as brown dwarf stars, whereas those secondaries that stay in the 
disc - and particularly those that stay in the inner disc - tend to 
grow by accretion, and sometimes they even grow larger than the 
primary. 

Switching from the standard barotropic equation of state to 
our new more realistic treatment of the energy equation has sev- 
eral systematic effects: 

- the efficiency of star formation (// = J {M+} /M C0RE ) is re- 
duced significantly (by ~ 15%); 

- the number of protostars formed (AC) is greatly increased 
(by -40%); 

- a higher proportion of brown dwarf stars is formed; 

- the mean period of multiple systems is reduced (by a factor 

~3); 

- the orbital eccentricities of multiple systems tend to be 
higher; 

- the mass ratios of multiple systems tend to be higher (i.e. 
more nearly equal components). 

All these trends can be attributed to the fact that the 
barotropic equation of state assumes that all gas is at the cen- 
tre of a collapsing spherical 1 M Q protostar, and therefore it be- 
comes adiabatic at relatively low densities. In contrast, our new 
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more realistic treatment of the energy equation allows the gas in 
low-mass proto-fragments to remain approximately isothermal 
to relatively high densities. 

These simulations do not capture all the deterministic effects 
in star formation. In particular, we do not take account of radia- 
tive and mechanical feedback from stars, we do not include non- 
ideal magneto-hydrodynamics effects, and we invoke sink parti- 
cles (which must prejudice interactions between stars, in favour 
of elastic dynamical ejections, and against dissipation and merg- 
ers). We plan to explore these additional effects in subsequent 
papers. In particular, we anticipate that by including feedback 
we can reduce the efficiency of star formation from the very 
high levels produced here ( ~ 60%) to values more compatible 
with observation ( < 30%; Alves et al. 2007; Nutter & Ward- 
Thompson 2007; Goodwin et al. 2008). 
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Appendix A: Smoothed mass distributions 

The smoothed mass distributions on Fig.|5]are given by a sum of 
Gaussians, 



P"** = I {(2^ exp 
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d/j 
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where fi = log 1() (M) and yi\ = log 10 (Mi). The standard deviation 
cr, is evaluated by adding - in quadrature - the mean separation 
between all masses across the entire mass spectrum (this is the 
first term on the righthand side of Eqn. IA.21 >. and the mean sep- 
aration between the five nearest masses (this is the second term 
on the righthand side of Eqn. |A.2t . Thus cr, combines a global 
and a local contribution to the smoothing. This smoothing is es- 
sentially ad hoc, and is designed purely to enable us to extract 
the large-scale features of the mass distribution, which are oth- 
erwise lost in the rather noisy histograms. In particular we are 
interested in the apparent bi-modality in the mass distributions 
obtained with the new treatment of the energy equation. 



